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A fast, fully Implicit approximate factorization algorithm designed to solve the conservative, transonic, full- 
potential equation in either two or three dimensions is described. The algorithm uses an upwind bias of the 
density coefficient for stability in supersonic regions. This provides an effective upwind difference of the 
streamwise terms for any orientation of the velocity vector (i.e., rotated differencing), thereby greatly enhancing 
the reliability of the present algorithm. A numerical transformation is used to establish an arbitrary body-fitted, 
finite-difference mesh. Computed results for both airfoils and simplified wings demonstrate substantial im- 
provement in convergence speed for the new algorithm relative to standard successive-line over-relaxation 
algorithms. 


Introduction 

A N implicit approximate factorization algorithm (AF2) 
for solving the low-frequency (unsteady) transonic small- 
disturbance equation in two and three dimensions was 
presented in Ref. 1. This algorithm has been subsequently 
applied to the solution of the transonic small-disturbance 
potential equation 2 and the conservative full-potential 
equation 3 * 5 for steady flows in two space dimensions. For 
both steady formulations, significant improvement in con- 
vergence speed has been obtained relative to the standard 
transonic relaxation procedure, successive-line over- 
relaxation (SLOR). In the present study, the AF2 algorithm is 
applied to the conservative full-potential equation for steady 
three-dimensional transonic flowfields. 

Several general guidelines for the construction of implicit 
approximate factorization (AF) schemes can be formulated by 
considering the two-level iteration procedure 

NC n +o>L<t> n =0 (1) 

where C n is the correction ( 4 >' I + / — <t> n ), L<t> n is the residual, 
which is a measure of how well the finite-difference equation 
is satisfied by the mh level velocity-potential solution (</>"), 
and 03 is a relaxation parameter. The iteration scheme given by 
Eq. (1) can be considered as an iteration in pseudotime, where 
the n superscript indicates the time step level of the solution; 
i.e., ( ) n+i — ( ) n - At ( ),. The operator N determines the 
type of iterative procedure, and therefore determines the rate 
at which the solution procedure converges. In the AF ap- 
proach, N is chosen as a product of several factors, usually 
two factors for two-dimensional algorithms and three factors 
for three-dimensional algorithms. These factors are chosen so 
that their product is an approximation to L, only simple 
matrix operations are required, and the overall scheme is 
stable. 

Stability in the present full-potential formulation for 
supersonic regions of flow has been achieved by the addition 
of an artificial viscosity term similar to that introduced in 
Ref. 6 . However, in the present formulation, addition of the 
artificial viscosity term is achieved by using an upwind bias of 
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the density coefficient. This strategy is significant because it 
simplifies the technique for including an artificial viscosity 
term into the residual operator. Other studies 79 have used 
similar steady-state differencing procedures in a wide variety 
of problems to further substantiate this differencing 
procedure as being both reliable and flexible. 

Governing Equations 

The three-dimensional full-potential equation written in 
strong conservation-law form is given by 

(p<t>x )* + (p<f> y ) y + (p<f> z ) z = 0 ( 2 a) 

[ 7 — 7 V/y-i 

J (2b) 

The density (p) and velocity components and <t> z ) are 

nondimensionalized by the stagnation density (p s ) and the 
critical sound speed (a*), respectively; x,y, and z are Cartesian 
coordinates in the streamwise, spanwise, and vertical 
directions, respectively, and 7 is the ratio of specific heats. 
The two-dimensional conservation-law form of the full- 
potential equation is simply obtained by dropping all y- 
derivative terms from Eq. (2). 

Equation (2) expresses mass conservation for flows that are 
steady, isentropic, and irrotational. The corresponding shock- 
jump conditions are valid approximations to the Rankine- 
Hugoniot relations for many transonic flow applications. A 
comparison of isentropic and Rankine-Hugoniot shock polars 
is given in Ref. 10. 

Equation (2) is transformed from the physical domain 
(Cartesian coordinates) to the computational domain by using 
a general independent variable transformation. This trans- 
formation, indicated by (see Fig. 1) 

Z = ${x,y,z) ti = y(x,y,z) t=t(x,y,z) (3) 

maintains the strong conservation-law form of Eq. (2). M The 
full-potential equation written in the computational domain 
coordinate system) is given by 

(tKtKt),- 

[ 7 — 7 1 1/y-1 

+ f )J (4b) 
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where 


U=A } <t>t +A 4 <f> rf +A s <t> ( A j -i 2 x + + f^ 

+A 2 <f> Tf +A 6 <j>{ A 4 = $ x ri x +Z y T) y +Z z ri z 
W-A^ k +A 6 $ ri +>4j4> f A s =f jr f r + {,{> + £zf z 


>!,={'+$'+{* 

^ =7 ?i4-^+T?^ 


^<5 =Vxfx + + 1 7 z f* 


(5a) 


and 


^jrVy fr £> *?z £r 4" %zVx £y %zVyi x %ytfx£z ZxVzty (^) 

U, V, and W are contravariant velocity components along the 
{, 17 , and f directions, respectively, A r A 6 are metric quan- 
tities, and J is the Jacobian of the transformation. To 
evaluate the expressions of Eq. (5), the following metric 
identities are necessary: 

^ x =J(y^-y ( z lt ) nx=J(y { z i -y i z t ) i x =J(y i z„-y n z t ) 

^=y(jf r z,-jf,z f ) r) y =J(X ( Z{-XtZ ( ) $ y =J(x y z i -x l z 1 ) 

tz =J(x 1 ,y t ) Vz =J(x i y i -x i y { ) = J(x t y v -x„y ( ) 


( 6 ) 

The two-dimensional form of the full-potential equation 
written in the computational domain ({-f coordinates) is 
obtained by dropping all y and jj terms in Eqs. (4-6); i.e., all y 
and tj derivatives as well as all derivatives of y and rj are set 
equal to zero. An exception to this is that y n and r\ y must be 
set equal to one. 


Several significant advantages are offered by this very 
general form. The main advantage is that boundaries 
associated with the physical domain are transformed to 
boundaries of the computational domain. This aspect is 
illustrated in Fig. 1 where the physical and computational 
domains for a typical transformation are shown. The com- 
putational coordinates {,rj, and f are in the wraparound, 
spanwise, and radial-like directions, respectively. The inner 
wing boundary transforms to f = f mtx , and the outer physical 
boundary transforms to f=f m j n . Note that no restrictions 
have been placed on the shape of the outer boundary. Ar- 
bitrarily shaped outer boundaries, including wind-tunnel 
walls, may be used. The symmetry-plane boundary trans- 
forms to v — Vmin* and the wing-tip boundary transforms to 
T; = 7 / max . The last two sides of the computational domain are 
formed from the upper and lower cuts along the vortex sheet. 

In the present study, the generality of Eqs. (4) and (5) is 
reduced somewhat by one simplification, namely, all 
ij = constant surfaces must coincide with y = constant planes. 
This can be expressed mathematically by 

y t =0 y ( =0 (7) 

The metric quantities of Eqs. (5) and ( 6 ) reduce to the 

following: 

Al=tl+t 2 y+£l A 2 =1) 2 y 

a + + A 4 =i y n y 

A, ={,{-, A 6 =r yVy ( 8 a) 

J= (fcrii -tx?x)v,=t/(x ( z { -x { z ( )y , ( 8b ) 


and 




BOUNDARY (r? = t? mjn ) 

ig. 1 Schematic of general <x t .y,z)~({,i|,f) transformation: a) 
lysical domain; b) computational domain. 


{,=/(*<.*„ -X n Z t ), t y= J (X 1J Zi “*{*,) 

S z = -Jx c y v , 

v x = o Vy-Uy^ v z =o ( 9 ) 

Use of Eq. (7) permits the simplified treatment of wing 
geometries and does not affect the generality of the three- 
dimensional spatial-differencing scheme or the fully implicit 
AF2 iteration scheme. 

Grid Generation 

The grid-generation scheme used in the present three- 
dimensional formulation is a simple extension of the two- 
dimensional scheme presented in Ref. 4. The finite-difference 
mesh for each spanwise plane (rj = constant plane) is generated 
using the standard two-dimensional algorithm. This requires 
solution of the following two Laplace equations in each 
spanwise plane: 

ixx+f zz=0 txx+tzz=0 < 10 > 

These equations are transformed to (and solved in) the 
computational domain; that is, { and f are the independent 
variables, and x and z are the dependent variables. A fast 
approximate factorization relaxation algorithm is used to 
solve the resulting transformed equations . 4 This establishes 
values for x and z in each spanwise plane. Coordinate values 
in the spanwise direction {y values) are established by a simple 
stretching formula (usually Ay = const). 

Next, the effects of wing sweep are built into the three- 
dimensional mesh by applying a simple translation to all x 
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values. This is accomplished by 

x ^ swept = x\ unswept "F {J 7)AytanA^ (11) 

where X, is the sweep angle at theyth mesh node in the span 
direction. Given the values of x, y, and z at all mesh points, 
the values of the quantities x t , x^ , etc., are easily computed 
using standard fourth-order-accurate finite-difference for- 
mulas. (It has been found that using fourth-order-accurate 
formulas, instead of second-order-accurate formulas, im- 
proves the accuracy of the resulting solutions.) Then, with 
Eqs. (8) and (9), the quantities £*,£,, and J are 
computed. 


where v is an artificial viscosity coefficient to be defined 
subsequently. The complete finite-difference approximation 
to Eq. (4a), including the addition of an upwind-differenced 
artificial viscosity term, is given by 



Full-Potential Equation Algorithm 
Spatial Differencing 

A second-order-accurate finite-difference approximation to 
Eq. (4a) is given by The , p„ , and derivatives of Eq. (16) have been evaluated 

with a backward (forward) difference when C/ l+ y jJtk , Vjj+x,ki 
and Wij'k+u are positve (negative), respectively. The r, s , and 
(12) t i n di ces control the difference direction on the density 
derivatives and are defined by 


s.(t) +J ’ ( 7 ) ,^( p t) 


y (pij,k+ Vi —PiJ,k + t+ l A (1^) 


where the i, j, and k subscripts indicate position in the £, 17, 
and f finite-difference mesh. The operators 

£*( ), S’, ( ), and £ f ( ) (13) 

are first-order accurate backward-difference operators in the 
£, rj, and f directions, respectively. The density calculation is 
performed in a straightforward manner by using Eq. (4b). 
Values of density are computed and stored at half-points in 
the finite-difference mesh; i.e., /+ l /i,j,k, in a manner similar 
to that suggested by Ref. 12. Values of <t> ^ , and <f>^ 

required for computing the density at /+ l AJ,k are given by 

i+t/jj'k = l// * (^i + lj+l.k ~ ^i+1 J-l,k +$ij+l.k ~4>iJ-l,k ) 
^ti+ttjk ~ (^i+U.k + l +/*>,*-/ *F $ij.k+l ^ij,k- 1 ) 

(14) 

The contravariant velocity components, U t+{/iJkt Vij+y ltk i 
and used both in Eq. (12) as well as in the density 

calculation, are computed by standard second-order accurate 
finite-difference formulas, an example of which is given by 


r— ±1 

when 

O, + ‘/ }Jt k s o 


s= ± 1 

when 

Kj+ Vi,k ^ 0 


/= ±1 

when 

Wjj'k+v, 

(18) 


This maintains an upwind influence in the differencing 
scheme for supersonic regions anywhere in the finite- 
difference mesh for any orientation of the velocity vector. 
Thus, use of the differencing scheme given by Eqs. (17) and 
(18) closely approximates the effects of a rotated differencing 
scheme. 5 This aspect greatly contributes to the stability and 
reliability of the present algorithm for many difficult test 
cases. 

The scheme given by Eqs. (17) and (18) is centrally dif- 
ferenced, hence, second-order accurate in subsonic regions. In 
supersonic regions, the differencing is a combination of the 
second-order accurate central differencing used in subsonic 
regions and the first-order accurate upwind differencing 
resulting from the addition of artificial viscosity. As the flow 
becomes increasingly supersonic, the scheme is increasingly 
retarded in the upwind direction. 

As shown in Refs. 3 and 4, Eq. (17) can be written in a 
simplified form given by 


Ui+HJ.k + ~4>U.k ) + l/4A 4 t +y t j,k (♦/+/./+/.* 

k -<t>ij-i,k) + ^^Sj+'AJ'k (^i+U.k + i 

* 

Calculation of the density at half-points in the finite- 
difference mesh produces better shock-wave resolution, first 
verified in two dimensions by South. 13 This is because the 
computational module in the streamwise £ direction extends 
over fewer grid points. Special formulas replace Eqs. (14) and 
(15) at boundaries and will be discussed in a subsequent 
section. 

Equation (12) is a suitable finite-difference scheme for 
subsonic flow regions. However, for supersonic regions, a 
properly chosen artificial viscosity term which usually 
provides an upwind bias to the differencing scheme, must be 
introduced to maintain stable convergence. 6 With the present 
formulation, stable supersonic regions can be maintained by 
adding an artificial viscosity term of the form 



The addition of the artificial viscosity given by Eq. (16) is thus 
equivalent to retarding the density in the original centrally 
differenced scheme [Eq. (12)]. Artificial viscosity is not added 
explicitly as in Ref. 6. The present approach is significant 
because it simplifies the technique for including an artificial 
viscosity term into the residual operator. 

The artificial viscosity coefficient, v , strongly affects the 
stability of the present scheme and is computed as follows: 


/ ll/l N 

, / IKI 

\ / 1 W\ \ 


„ r 


1 — Aij [vp v — 
{ x J 

j *-> 

£, 

< 

i 

(16) 

v i+ ViJ.k ~ j 


max[ (Mjj k —1)0,0] for U u l/iJtk > 0 
max [(M ?+ jj tk - 1)C,0] for U t +, /lJik <0 


( 20 ) 
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The parameter C is a user-specified constant and is usually set 
between 1.0 and 2.0. Use of larger values of C increases the 
amount of artificial viscosity and, therefore, the amount of 
upwinding in the difference scheme. An additional constraint 
is which improves the stability and in some cases, im- 
proves the convergence rate. Expressions of v at ij + l A,k and 
ij,k + l A are required and are written similarly to ^ J k . 


AF2 Iteration Scheme 

The AF2 fully implicit approximate factorization scheme 14 
is extended in the present study to the three-dimensional full- 
potential equation. This new scheme can be expressed by 
choosing the Noperator of Eq. (1) as follows: 

‘wcfc* =-[(«- MA ) (»* - *- MA ) 

-a£? , i4*](a + r r )CJ w (21a) 

where 


A,.(^y a,. 

\ J 'i-'AJ.k 



n 

ij- */j,k 


A k 



n 

ij.k- 


>/} 


( 21 b) 


The p, p, and p coefficients are defined by Eqs. (19b-d), a is a 
free parameter defined subsequently, and the operator Ef l is 
a shift operator given by 

>w* = ( )u.*+i ( 22 > 

Note that one form of the two-dimensional AF2 iteration 
scheme is obtained from Eq. (21a) by simply setting the 77 
difference equal to zero . 4 

Inclusion of the retarded density coefficients, p, p, and p, in 
the N operator is not necessary for stability. This situation 
was first discussed in Ref. 7 for an implicit ADI-type iteration 
scheme. Recent testing of the present AF2 iteration scheme in 
two dimensions, with the retarded density coefficients in the 
N operator replaced by just the density, yielded stable results 
with essentially no reduction in the convergence rate. Ad- 
ditional tests with the density coefficients in the N operator 
completely removed (i.e., replaced by one), also yielded stable 
results but slowed convergence by a factor of 2 to 3. 
Therefore, for the present formulation, existence of at least 
the density in the N operator, but not necessarily the up- 
winded or retarded density coefficients, is very important for 
achieving fast convergence but not always necessary for 
stability. 

Implementation of the AF2 scheme is achieved by writing it 
in a three-step form given by 


* Step 1 : 



T K A A)s?j= auL 'i > iM +aA k+lf?J.k+l 
A k 7 


(23) 


step 1 , the g array is obtained by solving a tridiagonal matrix 
equation for each £ = constant line in the /rth plane. In step 2, 
the / array is obtained from g by solving a tridiagonal matrix 
equation for each 77 = constant line, again for just the Arth 
plane. Next, step 1 is used to obtain the g array for the k a 1 
plane, and then step 2 is used to obtain the/ array for the k a 1 
plane, etc. This process continues until all values of / in the 
three-dimensional mesh are established. Then, by using step 
3, the correction array is obtained from the /array by solving 
a simple bidiagonal matrix equation for each f = constant line 
in the entire finite-difference mesh. The nature of this AF2 
factorization places a sweep-direction restriction on the step 
1-2 combination and on step 3. The step 1-2 combination must 
be swept in the direction of the decreasing k subscript; that is, 
from the wing boundary toward the outer boundary (see Fig. 
1). The step 3 sweep must proceed in just the opposite 
direction; that is, from the outer boundary toward the wing. 
There are no sweep-direction limitations placed on any of the 
three sweeps due to flow direction. 

Initiation of step 1 at the wing boundary ( k=NK ) requires 
knowledge of / at NK+ 1, which is generally unobtainable. A 
simple solution is to set / at NK + 1 equal to zero. Because the 
present iteration scheme is written in the correction form, / 
must approach zero as the solution converges. This boundary 
condition is therefore consistent with the steady-state 
solution. Other approaches using extrapolations of old known 
values of/(e.g.,/ 5 £ ; and could perhaps pro- 
vide faster convergence, but have not yet been investigated. A 
similar boundary condition is required for g at = 0=0 

and T 7 = i?max( i=NJ) and i s implemented by imposing 
(#! ,)/,/ = (£*) j.w = 

Temporal Damping 

For the AF2 factorization, the N operator must be written 
so that either the £-, 77-, or f-difference approximation to the 
full-potential equation is split between two factors. This 
construction generates either a 0 ( ,-, or 0 ^-type term, 
and if it is properly upwind differenced , 14 provides time- 
dependent dissipation to the convergence process. When a 
particular coordinate direction is split (e.g., the £ direction), 
the resulting <f> k( difference direction is fixed by the con- 
struction of the AF2 algorithm; that is, the <i> v term is either 
backward or forward differenced over the entire mesh. Due to 
the wraparound £ coordinate, a backward- (forward-) dif- 
ferenced <A f/ is upwind (downwind) differenced below the 
wing and downwind (upwind) differenced above. Therefore, 
a problem with <f> it arises either above or below the wing. 
Following the two-dimensional algorithm development , 4 the 
f-difference approximation is split between two factors. This 
allows control over the other more important coordinate 
directions (£ and 17) because the and <j> n , terms are added to 
the iteration scheme explicitly and are not part of the fac- 
torization construction. The <£„, and <t>^ { terms are included by 
adding 

=F|8, IK w *I^ and (26) 


Step 2: 

(24) 

\ a 7 

Step 3: 

(<x + f { )C? Jik =f?j, k (25) 

Here, is a relaxation factor equal to 1.8 for all cases 
presented; g'-j is an intermediate result stored at each grid 
point in a given k plane, i.e., g requires only a two- 
dimensional array of storage; and f” Jk is an intermediate 
result stored at each point in the finite-difference mesh. In 


inside the brackets of the first and second sweeps, Eqs. (23) 
and (24), respectively. The parameter is determined as 
follows: 


0* if 

0 if 


f M i+Uk >l 

r M i+IJtk < i 

La 


upper surface 
lower surface 

upper surface 
lower surface 


(27) 


where &^ M>i is a user-specified constant which can be ad- 
justed as needed. The parameter is a user-specified con- 
stant fixed over the entire mesh. For all cases presented in the 
present study, = 0. The double arrow notation on the £- and 
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tj - difference operators indicates that the difference is always 
upwind. For the y direction, a backward difference is used 
when the ij contravariant velocity component (V,j t k ) is 
positive, and a forward difference is used when V \s negative. 
The sign is chosen in each case so that the addition of <f> n{ and 
<f> it increases the magnitude of the first- and second-sweep 
diagonal coefficients, respectively. 

Of course, with this N operator construction, the 4> fr - type 
term is properly differenced in only the aft half of the mesh. 
In the two-dimensional algorithm presented in Ref. 4, the 
type term was properly differenced in the forward half of the 
mesh. Both versions seem to yield equivalent results. 
However, most experimentation has been conducted on the 
Ref. 4 version. Although adverse effects for cases with 
supersonic flow at the trailing edge may be anticipated with 
this form, none has been experienced. In fact, cases with 
freestream Mach numbers near unity have been computed (see 
the two-dimensional results section) in which the trailing edge 
is entirely imbedded in supersonic flow with no adverse ef- 
fects. 

The quantity a appearing in Eq. (21a) can be considered as 
At~ l . This direct analogy to time provides one strategy for 
obtaining fast convergence, namely, advance time as fast as 
possible with large time steps (i.e., small values of a). As 
pointed out in Refs. 2 and 3, this is effective for attacking the 
low-frequency errors but not the high-frequency errors. The 
best overall approach is to use an a sequence containing 
several values of a. The small values are particularly effective 
for reducing the low-frequency errors, and the large values are 
particularly effective for reducing the high-frequency errors. 
The a sequence given in Ref. 3 has been used for both two- 
and three-dimensional cases presented herein. The endpoints 
generally used in two dimensions are a L =0.07 and ot H = 1.5, 
and for three dimensions a L = 0.4 and a H = 4.0. 

Boundary Conditions 

The wing surface boundary condition is that of flow 
tangency (i.e., no flow through the wing surface), and 
requires the f contravariant velocity component at the wing 
surface be zero (W=0). This boundary condition is im- 
plemented by applying 



(28) 


where k=NK is the wing surface. In other expressions, where 
is required at the wing surface [Eqs. (14) and (15)], the 
W-0 boundary condition is used again to obtain 





(29) 


Thus, a value of 0 r at the wing surface can be obtained 
without using a one-sided difference on <£. 

In the present study, a special wing geometry has been 
chosen to evaluate the new three-dimensional AF2 algorithm, 
namely, flow past an arbitrary wing mounted between parallel 
walls. The purpose of this model problem is to simulate the 
flow past a wing in a wind tunnel. The parallel sidewalls are 
treated with the same tangency boundary condition used for 
the wing surface (i.e., K=0). 

Computed Results 

The implicit algorithm presented in the previous section has 
been coded into a transonic airfoil analysis computer code 
(TAIR) and a transonic wing analysis computer code 
(TWING). (For details of the TAIR code, see Ref. 4.) Each of 
these codes is evaluated in this section by presenting a range of 
computed examples. The two-dimensional results computed 
from TAIR were all computed in the default mode. This 
simply means that only three inputs are allowed to be changed 


from case to case: freestream Mach number, angle of attack, 
and airfoil coordinates. Other parameters including all 
relaxation factors, acceleration parameters, and temporal 
damping coefficients are either held fixed or are adjusted 
automatically by internal computer code logic. This feature 
greatly simplifies operation of TAIR and, at the same time, 
improves the reliability, especially for inexperienced users. In 
the default mode, the timelike dissipation coefficient 0 k and 
the artificial viscosity coefficient C are chosen very con- 
servatively =5, C = 2). Then as the solution progresses, 
both the number of supersonic points (NSP) and the cir- 
culation T are monitored. If NSP and T grow very rapidly 
indicating a difficult case, 0% and C are increased. Con- 
versely, if NSP and T grow slowly, 0 k and C are decreased. 
This philosophy generally produces converged results for 
most airfoil solutions ranging from subcritical cases to dif- 
ficult strong shock cases; however, it reduces the convergence 
speed below optimum by about 10-50^0, depending on the 
particular case. 

Two-dimensional comparisons with the GRUMFOIL 
computer code 15 are considered next. GRUMFOIL is similar 
to TAIR in that both codes solve the conservative full- 
potential equation, but different in that TAIR uses the AF2 
iteration scheme and GRUMFOIL uses a hybrid direct- 
solver/SLOR iteration scheme. 6 This hybrid iteration scheme 
is composed of one direct-solver iteration, which is very ef- 
fective for reducing low-frequency errors, but is unstable for 
supersonic regions, followed by several (10 is the default) 
SLOR iterations. The purpose of the SLOR iterations is to 
smooth high-frequency errors generated by the direct-solver 
step in regions of supersonic flow. (It should be pointed out 
that stable operation from the direct-solver iteration scheme, 
without benefit of the SLOR scheme, can be obtained for 
solutions in which the supersonic region is small, providing 
the level of artificial viscosity is increased above the standard 
amount. This aspect has not been investigated in the present 
study.) The boundary-layer option, which is available in 
GRUMFOIL, has not been used in any of the results 
presented herein. All GRUMFOIL results are computed on a 
148x32 mesh, while the TAIR results are computed on a 
149x32 mesh. 

The first two-dimensional case considered here is for the 
NACA 0012 airfoil at A/* =0.75 and a = 1 deg. The present 
surface pressure coefficient distribution for this calculation is 
compared in Fig. 2 with the GRUMFOIL computer code 
result. The agreement is excellent everywhere except at the 
shock wave where the disagreement in position is about one 
mesh-cell width. This error causes a 2Vi°7o error in lift 
coefficient. 

Convergence history curves for this calculation are 
presented in Fig. 3. As discussed in Ref. 3, use of the residual 
operator for comparing different iteration scheme con- 
vergence histories can produce quite misleading results. This 
is because the residual operator is heavily biased toward the 
high-frequency end of the error spectrum, and therefore, does 
not produce a fair comparison for algorithms which treat low- 
frequency errors differently. A much more suitable means of 
comparing convergence histories is to use error. The rms error 
in the airfoil surface pressure is used in the present study and 
is defined by 


/r W -i , Nl ^ 

E? ms =[[ i {c;-c Pi )’\l£cl t ] (30) 

where C n p . is the surface pressure coefficient at the ith grid 
point ani the nth iteration, C Pi is the surface pressure 
coefficient at the ith grid point taken from the converged 
solution, and NI is the total number of surface grid points. 

The three curves shown in Fig. 3 correspond to the 
following iteration schemes: 1) AF2, 2) hybrid, and 3) SLOR. 
Each convergence-history curve is constructed by plotting 
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Fig. 2 Pressure coefficient comparison (NACA 0012 airfoil, 
A/« = 0.75, a=l deg). 



Fig. 3 Two-dimensional convergence histories (NACA 0012 airfoil, 
M m =0.75, « = 1 deg). 


E rms vs CPU time (Ames CDC 7600 computer). The SLOR 
scheme is simply the hybrid scheme without benefit of the 
direct-solver step, and has been approximately optimized by a 
trial-and-error adjustment of the relaxation parameters. The 
hybrid cases have been computed with default values for all 
relaxation parameters. Setup times, that is, the CPU time 
required for grid generation, initialization, and coarse- and 
medium-mesh calculations, are included in each convergence- 
history curve. The AF2 curve includes 6 s for grid generation 
and initialization. The hybrid and SLOR curves both use 
coarse-medium-fine mesh sequences. Converged results from 
the coarse mesh are interpolated onto the medium mesh, and 
then from the medium mesh onto the fine mesh, thus 
providing a good initial guess for the fine-mesh calculation. 
The setup times for these cases are 9 s for the hybrid scheme 
and 43 s for the SLOR scheme. For this calculation, AF2 is 
3.5 times faster than the hybrid scheme and 7 times faster than 
SLOR. 



deg). 


Solutions with Af* Approaching Unity 

As the freestream Mach number approaches one, in- 
teresting airfoil shock wave patterns develop. For instance, 
the shock wave pattern about an NACA 0012 airfoil at 
^^=0.95 and a = 4 deg is shown by the Mach number 
contours in Fig. 4. A so-called “fishtail” shock system is 
formed. Relatively weak supersonic-to-supersonic oblique 
shocks emanate from the trailing edge and merge with a 
normal shock downstream of the airfoil. The oblique shock 
emanating from the trailing edge upper surface has been 
strengthened by the addition of circulation, while the oblique 
shock emanating from the trailing edge lower surface has been 
weakened and is almost nonexistent. The normal shock above 
the airfoil is much stronger than the normal shock below the 
airfoil. The triangular region between the oblique and normal 
shocks has a nearly constant Mach number which is ap- 
proximately 1.1. This shock wave pattern is characteristic of 
solutions with freestream Mach numbers near unity and has 
been observed experimentally as well as computationally. It is 
generally considered to be the correct qualitative solution. For 
instance, a fishtail shock solution for a 10°7o circular arc 
airfoil at A/* = 0.98 and 0 deg angle of attack was presented in 
Ref. 7. This calculation was a solution to the conservative 
full-potential equation using a Cartesian mesh and small- 
disturbance boundary conditions. Because the flow in this 
solution is essentially aligned with the finite-difference mesh, 
“rotated differencing” is not necessary. 3 However, with the 
present wraparound coordinate system, rotated differencing 
is essential for maintaining stability of difficult cases such as 
the present fishtail shock solution. 5 

As pointed out in Ref. 14, nonconservative full-potential 
equation solutions with freestream Mach numbers ap- 
proaching unity are characterized by strong oblique shock 
waves at the trailing edge followed by subsonic flow. The 
fishtail shock structure for these cases is not predicted. 
Conservative vs nonconservative differencing was the subject 
of discussion in Ref. 16 where similar differences for the 
transonic small-disturbance potential equation were reported. 
It is generally understood that these differences are the result 
of effective mass creation at shock waves for the non- 
conservative differencing schemes. Therefore, to obtain the 
proper mass balance and the correspondingly correct solution 
conservative differencing is required. 

Convergence history curves for the last calculation, in- 
cluding E rms , I/? I max [defined by Eq. (19a)], NSP, and Q 
convergence histories, are presented in Fig. 5. Convergence is 
achieved in this case in approximately 20 s of CPU time (116 
iterations), as indicated by constant values of both NSP and 
C L . At this point, l/?l max has dropped only slightly, while 
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Fig. 5 Two-dimensional convergence histories (NACA 0012 airfoil, 
Af„ = 0.95, a = 4 deg). 



x/c 

Fig. 6 Two- and three-dimensional pressure coefficient comparison 
(NACA 00XX airfoil). 


£ rrm has dropped by over 2 Vi orders of magnitude. This 
discrepancy is partially due to the fact that the airfoil surface 
• solution, which is effectively monitored by £ rms , converges 
before the downstream fishtail shock wave pattern. Another 
cause is that during the initial phase of convergence in which 
the residual does not drop, the position of the shock sonic line 
is rapidly being adjusted. This excites high-frequency errors, 
and therefore keeps the residual artificially high even though 
E is being reduced. Establishment of convergence for such 
a small reduction in the maximum residual is a characteristic 
behavior of many strong shock calculations using the present 
algorithm. Use of a three order-of-magnitude reduction in 
\R\ as a convergence criterion for the present case would 
require £ rms to be reduced by five orders of magnitude, which 
represents a factor of two more iterations than necessary. 

This is the first time in transonic flow computations that 
calculations such as the last one have been computed using the 
conservative full-potential equation with an exact airfoil 
mapping. Without the newly developed rotated difference 
scheme, these calculations would have been unstable. The 



Fig. 7 Shock wave sonic line position relative to the wing planform 
(NACA 0015 wing, Af* = 0.S6, X = 30 deg, /R = 1.9, a = 0 deg). 

rapid convergence of this difficult case demonstrates the 
reliability and efficiency of the present transonic flow solution 
procedure. 

Three-Dimensional Solutions 

Results from the transonic wing analysis code (TWING) are 
presented in this section. All calculations have been computed 
with the density coefficients upwinded in only the £ and rj 
directions. For the cases presented herein, this was sufficient; 
but for cases with stronger shocks at the trailing edge, density 
upwinding along the f direction would probably be required. 
All results have been computed with ranging from 0.1 to 
0.5. The larger values of were required for the larger aspect 
ratio cases. 

To evaluate the three-dimensional code, several infinite- 
aspect-ratio results have been compared with two- 
dimensional results using the concept of simple sweep theory. 
For three-dimensional infinite-aspect-ratio calculations, the 
solution in the wing-normal plane is the same as a two- 
dimensional solution with appropriate scaling; that is, 
M m 2D = cos\-A/„ fJD , C Pt3D =cos 2 \-C Pf2D , and 
T(x) 2D = T(x) 3D /c osX, where T(x) is the airfoil thickness 
distribution. Subcritical results using this concept have been 
compared and are in excellent agreement. A supercritical 
comparison is shown in Fig. 6. The three-dimensional 
calculation was computed with an aspect ratio (/R) of 9.5 and 
thus simulates an infinite-aspect-ratio condition at center 
span. Other conditions for this calculation were NACA 0015 
wing, Af* =0.86, X = 25 deg, and a = 0 deg. The two results 
are in close agreement but are not identical. The reason for 
the slight disagreement is numerical and is due to the in- 
creased level of artificial viscosity arising from larger values 
of the local Mach number in the three-dimensional case. If the 
artificial viscosity parameter (v) in the two-dimensional case 
is increased by adding sin^ \Af £ t3D to the M 2 term of Eq. (20), 
the two numerical solutions become identical. 

Figures 7 and 8 display the results of a wing calculation 
with the following characteristics: NACA 001 5 wing, yR = 1 .9, 
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Fig. 8 Three-dimensional pressure coefficient distribution (NACA 
0015 wing, = 0.86, A = 30 deg, /R = 1.9, a = 0 deg). 



Fig. 9 Center span pressure coefficient evolution with iteration 
number (NACA 0015 wing, /R=1.9, a = 0 deg, =0.86, A = 30 
deg). 


A/«, =0.86, A = 30 deg, and a = 0 deg. The shock wave sonic- 
line position is plotted to scale on the wing planform in Fig. 7. 

*As expected, the shock wave approaches both sidewalls in a 
nearly perpendicular fashion. Between the sidewalls near 
center span, the shock is swept and is approximately parallel 
to the wing leading and trailing edges. The shock-strength 
gradient along the span is quite large as indicated by the 
surface C p distributions shown in Fig. 8. The maximum local 
Mach numbers in the root, center span, and tip planes are 
1.26, 1.33, and 1.78, respectively. A large part of this 
spanwise shock-strength gradient is caused by the tip 
sidewall/wing interaction, which is essentially the opposite of 
a three-dimensional relief effect. The existence of the tip 
sidewall constrains the streamlines to remain in the tip plane. 
The wing sweep induces a spanwise component of velocity 
which in effect squeezes the streamlines toward the tip plane. 
This increases the Mach number and therefore the shock 
strength in the tip plane region. 

Figure 9 displays the center span C p distribution evolution 
with iteration number (n) for the case just presented. The 
n = 0 solution corresponds to the initial condition, and the 
n- 146 solution corresponds to the solidly converged solution 


in which the maximum residual has been reduced by three 
orders of magnitude. In just 20 iterations, a good ap- 
proximation is established in which the shock position is off 
by only l°7o. By 60 iterations the difference in shock positions 
is less than \°7o. 

All TWING calculations presented herein have been 
computed using a mesh with 46,000 points (115x20x20). 
This places 58 points from leading to trailing edge on both the 
upper and lower surfaces for each spanwise plane. Each 
iteration requires about 4.2 s of CPU time on the Ames CDC 
7600 computer. For the case of Fig. 8, this equates to a total 
run time of about 5 min to achieve plottable accuracy. For 
subcritical cases, the run time is much less, amounting to just 
over a minute to achieve plottable accuracy. 

Conclusions 

A fast, implicit algorithm for solving the conservative full- 
potential equation in both two and three dimensions is 
presented. Stability in supersonic regions is maintained by 
using an upwind evaluation of the density coefficient along all 
coordinate directions. This provides an effective upwind 
difference of the streamwise terms for any orientation of the 
velocity vector (i.e., rotated differencing), and thereby greatly 
enhances the reliability of the present algorithm. Use of the 
newly developed rotated differencing scheme has been in- 
strumental in computing a number of difficult two- 
dimensional test cases including several cases with “fishtail” 
shock-wave patterns. This represents the first time such 
calculations have been computed using the conservative full- 
potential equation with an exact airfoil mapping. 

The present fully implicit AF2 algorithm has been com- 
pared with both the standard transonic-solution procedure, 
successive-line over-relaxation (SLOR) and a hybrid (direct- 
solver/SLOR) scheme. The surface C p distributions produced 
by these schemes are in good agreement. Based on CPU time, 
the rms error in the surface pressures is reduced from five to 
seven times faster by the AF2 algorithm relative to SLOR. 
The hybrid scheme displays a very wide range of convergence 
speeds, being very fast for subcritical cases but much slower 
for strong shock cases. The AF2 scheme displays a much more 
uniform convergence rate over a very wide range of cases. It is 
essentially comparable to the hybrid scheme for subcritical 
cases but as much as 3.5 times faster for cases with strong 
supersonic regions. 

Three-dimensional results from the newly developed 
transonic-wing-analysis code (TWING) are presented. In- 
finite-aspect-ratio results are in good agreement with standard 
two-dimensional results. Other calculations for swept wings 
mounted between parallel walls are presented in which a 
strong shock wave extends across the entire wing span in- 
dicating a high degree of reliability. Convergence histories 
indicate that the substantial improvement in convergence rate 
experienced in two-dimensional cases carries over to the newly 
implemented three-dimensional version of the fully implicit 
AF2 iteration scheme. 
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